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Abstract. An efficient procedure for tfie computation of the coefficients of 
Legendre expansions is here presented. We prove that the Legendre coefficients 
associated with a function f{x) can be represented as the Fourier coefficients of 
an Abel-type transform of f{x). The computation of TV Legendre coefficients 
can then be performed in 0{N log N) operations with a single Fast Fourier 
Transform of the Abel-type transform of f{x). 



1. Introduction 

The efficient computation of the coefficients of Legendre expansions is a very im- 
portant problem in numerical analysis and applied mathematics with a wide range 
of applications including, just to mention a few, approximation theory, solution 
of partial differential equations and quadratures. Recently its relevance emerged 
also in connection with the computation of spectra of highly oscillatory Fredholm 
integral operators, which play an important role in laser engineering [2]. 

The difficulty of the problem lies essentially in the fact that these coefficients are 
represented by integrals whose integrands oscillate rapidly for large values of the 
index of the polynomials. Using standard quadrature procedures for the calculation 
of N Legendre coefficients leads only to slow 0{N'^) algorithms (see, e.g., Ref. [5 ). 
More efficiently, in Ref. [T] (see also [T31II1]) the Legendre coefficients are obtained 
by a suitable transformation of the corresponding Chebyshev coefficients, which 
yields faster 0{N{\ogN)'^) algorithms. Recently this problem has been also clearly 
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discussed in a paper by A. Iserles |10| . in which an algorithm for the computation 
of the Legendre coefficients, which is certainly fast and brilliant, is presented. 

In this paper we present an alternative procedure. The basic idea of our method 
consists in exploiting the Dirichlet-Murphy integral representation of the Legendre 
polynomials. Next, we prove that the coefficients of the Legendre expansion of a 
function f{x) are connected with a subset of the Fourier coefficients (the ones with 
nonnegative index) of an Abel-type transform of f{x). 

The numerical implementation of the algorithm follows straightforwardly and is 
very efficient. The aforementioned Fourier coefficients (which represent the searched 
Legendre coefficients) can be computed in 0{N logA^) operations by a single Fast 
Fourier Transform after the evaluation of the Abel-type integral by means of stan- 
dard quadrature techniques. 

2. Connection of Legendre expansions to Fourier series 
The standard form of the Legendre expansion reads: 

oo 

(2.1) /(x) = ^c„P„(x) xe[-l,l], 

where Pn{x) are the Legendre polynomials, which can be defined by the generating 
function [6]: 

(2.2) ^P„(x)r = (l-2a;t + t2) \ 

n=0 

and the coefficients {c„}^q are given by: 

(2.3) Cn=(^n+^ j J{x)P^{x)dx (n^O). 

The conditions to be satisfied by f{x) to guarantee the uniform convergence of 
the series in (|2.ip are discussed in [5]. However, for our purpose of computing the 
Legendre coefficients c„ it is sufficient to assume that f{x) be summable in the 
interval [—1,1]. 

We can now state the following theorem. 

Theorem 1. The coefficients {an}^=0' defined as an = (2n+i) ' coincide with the 
Fourier coefficients (with 0) of an Abel-type transform of f(x), that is: 

(2.4) a„ = = r f{y) e™^ dy (n > 0), 

(271+1) 

where the 2-k -periodic function f[y) is defined by 

(2.5) Ky)^^e{y)e^i f ^^""^ dx (y £ M), 

27ri J cosy [2(a; - cosy)]2 

£{y) being the sign function. 

Proof. Plugging the Dirichlet-Murphy integral representation of the Legendre poly- 
nomials ilSj Ch. Ill, §5.4]: 

(2.6) P„(cosa;) = / -dy, 

Jx [2(cosa: — cosy)] 2 
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into equality p.3|) (after the change of variable x cos a;), we have: 

(2.7) 2'K\an— I da; /(cos a;) sin X / r dy- 

"'o Jx [2(cosa; — cosy)]^ 

Interchanging the order of integration in ()2.7p we have: 

27ria„= rdye'("+^)^ f f{cosx) — -Ax 

Jo Jo [2(cosa; - cos?/)]2 

+ dye'("+^)W /(cosx) — -dx. 

Jtt Jo [2(cos a; — cos j/)] ^ 

Next, if we make the change of variables: y -> y — 27r and x — s> —a:, the second 
integral on the r.h.s. of (|2.8I) becomes: 

dye'("+^>y / /(cosx) — rdx. 

-TT "'o [2(cosa; — cosy)] ^ 

Finally, we obtain: 

27ria„= / dye'("+5)?' f f{cosx) ^— r da; 

"'0 Jo [2(cosx — cosw)l 2 

(2-10) ,0 . 

dye'("+5)W /(cosx) — -dx, 

Jo [2(cosa; — cosy)] ^ 

which, after the change of variable cos a; — ?> x into the integrals on the r.h.s., yields: 



(2.11) a„= / /(y)e'"^dy (n ^ 0), 

J —IT 

with 7(y) given by (^3]) . □ 
It is easy to check from (|2.5p that /(y) satisfies the following symmetry relation: 

(2.12) /(y) = -e'^/(-y). 



This latter, along with formulae \2A\ and ()2.5p . allows us to write the Legendre 
coefficients Cn in the following form: 



(2.13) c„ = 3 (n+i) / 0(y)sin 

where 



dy. 



(2.14) Hy)^ f — — rd^- 

J cos y [2{x - cosy)] 2 

The numerical implementation of the algorithm first requires the computation of 
the Abel-type integral /(y) defined in (|2.5p (or, equivalently, of the function 0(y) 
in (j2.14p V The integrand presents a weak algebraic singularity at the end point of 
the domain of integration, which can be effectively handled by means of a proper 
nonlinear change of variable. This technique, along with the use of a standard 
quadrature formula (e.g., the Gauss-Legendre one), allows obtaining high accuracy 
with a small number of nodes [12] . 
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Table 1. True and computed Legendre coefficients c„ for the 
function f{x) = \x\^^^. 



n 


True Cn 


Computed c„ 


Error 








40ooooooooooonnoonnn 





400000000002681 87694 


2 


68- 


12 


2 


n 


UUUUUUUUUUUUUUUUUU'JU 


n 


66666666671 084839901 


4 


42 


11 


4 


-0 

\J 


0Q9307fiQ9'^07fiQ9307fiQ 


-0 

\J 


0Q2307fiQ24fi0294091 fiQ 


I 


53 


10 


6 




03Q91 'i686'?74'inQ80'^Q2 


yj 


03921 5686307690289^1 


'i 

O 


32 


11 


8 


-0 


021 Q78021 Q78021 Q7802 


-0 


021 978022261 84656509 


2 


84_ 


10 


10 





01411764705882352941 





01411764762633771833 


5 


68_ 


10 


12 


-0 


00985221674876847290 


-0 


00985221698441776129 


2 


36_ 


10 


14 





00727272727272727272 





00727272772068757959 


4 


48_ 


10 


16 


-0 


00559179869524697110 


-0 


00559179938291245667 


6 


88_ 


10 


18 





00443458980044345898 





00443458980116155730 


7 


18_ 


13 


20 


-0 


00360360360360360360 


-0 


00360360382001508440 


2 


16_ 


10 


22 





00298656047784967645 





00298656131882347585 


8 


41_ 


10 


24 


-0 


00251572327044025157 


-0 


00251572347582495201 


2 


05- 


10 


26 





00214822771213748657 





00214822836757538881 


6 


55_ 


10 


28 


-0 


00185586142901330034 


-0 


00185586305079012375 


1 


62_ 


09 


30 





00161943319838056680 





00161943394840411781 


7 


50- 


10 



Finally, formula p.4p makes it possible to take full advantage of the computa- 
tional efficiency of the Fast Fourier Transform both in terms of speed of computation 
and of accuracy [31 [U] ■ The calculation of the first N coefficients of the expansion 
can consequently be accomplished in 0{NlogN) operations. 

The algorithm described has been implemented in double precision arithmetics 
using the open source GNU Scientific Library (GSL) T, and its performance has 
been tested on a variety of functions. 

First, feasibility and accuracy of the algorithm have been verified by direct com- 
parison of the obtained numerical results with the true Legendre coefficients for the 
function /(x) = |a;p/^, whose Legendre coefficients are known to be |4i p. 78]: 







(2.15) c„ + 

(2n + l)a{a 



{a 



{a + l){a + 3)---{a + n + l) 



if n odd, 
if n = 0, 

otherwise. 



where a — 3/2. Values of the computed Legendre coefficients along with the abso- 
lute error are given in Table [1] 
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The increment of performances, with respect to the computation of the Legendre 
coefficients c„ by ordinary quadrature, has been verified in terms of speed of com- 
putation at (nearly) equahty of precision. All accuracies have been determined by 
comparing the results of the algorithm with the reference values of c„, computed 
with 20 significant figures by standard quadrature with Mathematica [iTj . For 
these tests we used various functions (many of them have been already used in 
previous works), including polynomials, exponential/hyperbolic functions, rational 
functions (e.g., f{x) — ^l'^^2 with 7 = constant). All the results have confirmed 
the enormous increase of computational speed (the expected improvement ratio be- 
ing proportional to TV/ log TV). Such an increase of performances will become even 
more crucial for the efficient evaluation of multivariate Legendre transform [2] and 
expansions in Gegenbauer (alias Ultraspherical Legendre) polynomials, which will 
be the subject of a forthcoming paper. 
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